{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4ff7325d6fa34f1fb74c89ff3ab76cb3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Plot(antialias=3, axes=['x', 'y', 'z'], background_color=16777215, camera=[4.5, 4.5, 4.5, 0.0, 0.0, 0.0, 1.0, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Solving Poisson equation using\n",
    "a mixed (two-field) formulation.\n",
    "\"\"\"\n",
    "# needs python3\n",
    "# https://fenicsproject.org/docs/dolfin/2018.1.0/python/demos/mixed-poisson\n",
    "from dolfin import *\n",
    "\n",
    "# Create mesh\n",
    "mesh = UnitSquareMesh(30, 30)\n",
    "\n",
    "# Define finite elements spaces and build mixed space\n",
    "BDM = FiniteElement(\"BDM\", mesh.ufl_cell(), 1)\n",
    "DG  = FiniteElement(\"DG\", mesh.ufl_cell(), 0)\n",
    "W   = FunctionSpace(mesh, BDM * DG)\n",
    "\n",
    "# Define trial and test functions\n",
    "(sigma, u) = TrialFunctions(W)\n",
    "(tau, v)   = TestFunctions(W)\n",
    "\n",
    "# Define source function\n",
    "f = Expression(\"10*exp(-(pow(x[0]-0.5, 2) + pow(x[1]-0.5, 2))/0.02)\", degree=2)\n",
    "\n",
    "# Define variational form\n",
    "a = (dot(sigma, tau) + div(tau) * u + div(sigma) * v) * dx\n",
    "L = -f * v * dx\n",
    "\n",
    "# Define function G such that G \\cdot n = g\n",
    "class BoundarySource(UserExpression):\n",
    "    def __init__(self, mesh, **kwargs):\n",
    "        self.mesh = mesh\n",
    "        super().__init__(**kwargs)\n",
    "    def eval_cell(self, values, x, ufc_cell):\n",
    "        cell = Cell(self.mesh, ufc_cell.index)\n",
    "        n = cell.normal(ufc_cell.local_facet)\n",
    "        g = sin(5 * x[0])\n",
    "        values[0] = g * n[0]\n",
    "        values[1] = g * n[1]\n",
    "    def value_shape(self):\n",
    "        return (2,)\n",
    "G = BoundarySource(mesh, degree=2)\n",
    "\n",
    "# Define essential boundary\n",
    "def boundary(x):\n",
    "    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS\n",
    "bc = DirichletBC(W.sub(0), G, boundary)\n",
    "\n",
    "# Compute solution\n",
    "w = Function(W)\n",
    "solve(a == L, w, bc)\n",
    "(sigma, u) = w.split()\n",
    "\n",
    "\n",
    "########################################################### vtkplotter\n",
    "from vtkplotter.dolfin import plot\n",
    "\n",
    "# Plot solution on mesh, and warp z-axis by the scalar value\n",
    "plot(u, warpZfactor=1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
